home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / zmul.c < prev   
Encoding:
C/C++ Source or Header  |  1996-09-07  |  29.5 KB  |  1,106 lines
  1. /*
  2.  * Copyright (c) 1994 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Faster than usual multiplying and squaring routines.
  7.  * The algorithm used is the reasonably simple one from Knuth, volume 2,
  8.  * section 4.3.3.  These recursive routines are of speed O(N^1.585)
  9.  * instead of O(N^2).  The usual multiplication and (almost usual) squaring
  10.  * algorithms are used for small numbers.  On a 386 with its compiler, the
  11.  * two algorithms are equal in speed at about 100 decimal digits.
  12.  */
  13.  
  14. #include "zmath.h"
  15.  
  16.  
  17. LEN _mul2_ = MUL_ALG2;        /* size of number to use multiply algorithm 2 */
  18. LEN _sq2_ = SQ_ALG2;        /* size of number to use square algorithm 2 */
  19.  
  20.  
  21. static HALF *tempbuf;        /* temporary buffer for multiply and square */
  22.  
  23. static LEN domul MATH_PROTO((HALF *v1, LEN size1, HALF *v2, LEN size2, HALF *ans));
  24. static LEN dosquare MATH_PROTO((HALF *vp, LEN size, HALF *ans));
  25.  
  26.  
  27. /*
  28.  * Multiply two numbers using the following formula recursively:
  29.  *    (A*S+B)*(C*S+D) = (S^2+S)*A*C + S*(A-B)*(D-C) + (S+1)*B*D
  30.  * where S is a power of 2^16, and so multiplies by it are shifts, and
  31.  * A,B,C,D are the left and right halfs of the numbers to be multiplied.
  32.  */
  33. void
  34. zmul(z1, z2, res)
  35.     ZVALUE z1, z2;        /* numbers to multiply */
  36.     ZVALUE *res;        /* result of multiplication */
  37. {
  38.     LEN len;        /* size of array */
  39.  
  40.     if (ziszero(z1) || ziszero(z2)) {
  41.         *res = _zero_;
  42.         return;
  43.     }
  44.     if (zisunit(z1)) {
  45.         zcopy(z2, res);
  46.         res->sign = (z1.sign != z2.sign);
  47.         return;
  48.     }
  49.     if (zisunit(z2)) {
  50.         zcopy(z1, res);
  51.         res->sign = (z1.sign != z2.sign);
  52.         return;
  53.     }
  54.  
  55.     /*
  56.      * Allocate a temporary buffer for the recursion levels to use.
  57.      * An array needs to be allocated large enough for all of the
  58.      * temporary results to fit in.  This size is about twice the size
  59.      * of the largest original number, since each recursion level uses
  60.      * the size of its given number, and whose size is 1/2 the size of
  61.      * the previous level.  The sum of the infinite series is 2.
  62.      * Add some extra words because of rounding when dividing by 2
  63.      * and also because of the extra word that each multiply needs.
  64.      */
  65.     len = z1.len;
  66.     if (len < z2.len)
  67.         len = z2.len;
  68.     len = 2 * len + 64;
  69.     tempbuf = zalloctemp(len);
  70.  
  71.     res->sign = (z1.sign != z2.sign);
  72.     res->v = alloc(z1.len + z2.len + 1);
  73.     res->len = domul(z1.v, z1.len, z2.v, z2.len, res->v);
  74. }
  75.  
  76.  
  77. /*
  78.  * Recursive routine to multiply two numbers by splitting them up into
  79.  * two numbers of half the size, and using the results of multiplying the
  80.  * subpieces.  The result is placed in the indicated array, which must be
  81.  * large enough for the result plus one extra word (size1 + size2 + 1).
  82.  * Returns the actual size of the result with leading zeroes stripped.
  83.  * This also uses a temporary array which must be twice as large as
  84.  * one more than the size of the number at the top level recursive call.
  85.  */
  86. static LEN
  87. domul(v1, size1, v2, size2, ans)
  88.     HALF *v1;        /* first number */
  89.     LEN size1;        /* size of first number */
  90.     HALF *v2;        /* second number */
  91.     LEN size2;        /* size of second number */
  92.     HALF *ans;        /* location for result */
  93. {
  94.     LEN shift;        /* amount numbers are shifted by */
  95.     LEN sizeA;        /* size of left half of first number */
  96.     LEN sizeB;        /* size of right half of first number */
  97.     LEN sizeC;        /* size of left half of second number */
  98.     LEN sizeD;        /* size of right half of second number */
  99.     LEN sizeAB;        /* size of subtraction of A and B */
  100.     LEN sizeDC;        /* size of subtraction of D and C */
  101.     LEN sizeABDC;        /* size of product of above two results */
  102.     LEN subsize;        /* size of difference of halfs */
  103.     LEN copysize;        /* size of number left to copy */
  104.     LEN sizetotal;        /* total size of product */
  105.     LEN len;        /* temporary length */
  106.     HALF *baseA;        /* base of left half of first number */
  107.     HALF *baseB;        /* base of right half of first number */
  108.     HALF *baseC;        /* base of left half of second number */
  109.     HALF *baseD;        /* base of right half of second number */
  110.     HALF *baseAB;        /* base of result of subtraction of A and B */
  111.     HALF *baseDC;        /* base of result of subtraction of D and C */
  112.     HALF *baseABDC;        /* base of product of above two results */
  113.     HALF *baseAC;        /* base of product of A and C */
  114.     HALF *baseBD;        /* base of product of B and D */
  115.     FULL carry;        /* carry digit for small multiplications */
  116.     FULL carryACBD;        /* carry from addition of A*C and B*D */
  117.     FULL digit;        /* single digit multiplying by */
  118.     HALF *temp;        /* base for temporary calculations */
  119.     BOOL neg;        /* whether imtermediate term is negative */
  120.     register HALF *hd, *h1=NULL, *h2=NULL;    /* for inner loops */
  121.     SIUNION sival;        /* for addition of digits */
  122.  
  123.     /*
  124.      * Trim the numbers of leading zeroes and initialize the
  125.      * estimated size of the result.
  126.      */
  127.     hd = &v1[size1 - 1];
  128.     while ((*hd == 0) && (size1 > 1)) {
  129.         hd--;
  130.         size1--;
  131.     }
  132.     hd = &v2[size2 - 1];
  133.     while ((*hd == 0) && (size2 > 1)) {
  134.         hd--;
  135.         size2--;
  136.     }
  137.     sizetotal = size1 + size2;
  138.  
  139.     /*
  140.      * First check for zero answer.
  141.      */
  142.     if (((size1 == 1) && (*v1 == 0)) || ((size2 == 1) && (*v2 == 0))) {
  143.         *ans = 0;
  144.         return 1;
  145.     }
  146.  
  147.     /*
  148.      * Exchange the two numbers if necessary to make the number of
  149.      * digits of the first number be greater than or equal to the
  150.      * second number.
  151.      */
  152.     if (size1 < size2) {
  153.         len = size1; size1 = size2; size2 = len;
  154.         hd = v1; v1 = v2; v2 = hd;
  155.     }
  156.  
  157.     /*
  158.      * If the smaller number has only a few digits, then calculate
  159.      * the result in the normal manner in order to avoid the overhead
  160.      * of the recursion for small numbers.  The number of digits where
  161.      * the algorithm changes is settable from 2 to maxint.
  162.      */
  163.     if (size2 < _mul2_) {
  164.         /*
  165.          * First clear the top part of the result, and then multiply
  166.          * by the lowest digit to get the first partial sum.  Later
  167.          * products will then add into this result.
  168.          */
  169.         hd = &ans[size1];
  170.         len = size2;
  171.         while (len--)
  172.             *hd++ = 0;
  173.  
  174.         digit = *v2++;
  175.         h1 = v1;
  176.         hd = ans;
  177.         carry = 0;
  178.         len = size1;
  179.         while (len >= 4) {    /* expand the loop some */
  180.             len -= 4;
  181.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  182.             *hd++ = sival.silow;
  183.             carry = sival.sihigh;
  184.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  185.             *hd++ = sival.silow;
  186.             carry = sival.sihigh;
  187.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  188.             *hd++ = sival.silow;
  189.             carry = sival.sihigh;
  190.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  191.             *hd++ = sival.silow;
  192.             carry = sival.sihigh;
  193.         }
  194.         while (len--) {
  195.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  196.             *hd++ = sival.silow;
  197.             carry = sival.sihigh;
  198.         }
  199.         *hd = (HALF)carry;
  200.  
  201.         /*
  202.          * Now multiply by the remaining digits of the second number,
  203.          * adding each product into the final result.
  204.          */
  205.         h2 = ans;
  206.         while (--size2 > 0) {
  207.             digit = *v2++;
  208.             h1 = v1;
  209.             hd = ++h2;
  210.             if (digit == 0)
  211.                 continue;
  212.             carry = 0;
  213.             len = size1;
  214.             while (len >= 4) {    /* expand the loop some */
  215.                 len -= 4;
  216.                 sival.ivalue = ((FULL) *h1++) * digit
  217.                     + ((FULL) *hd) + carry;
  218.                 *hd++ = sival.silow;
  219.                 carry = sival.sihigh;
  220.                 sival.ivalue = ((FULL) *h1++) * digit
  221.                     + ((FULL) *hd) + carry;
  222.                 *hd++ = sival.silow;
  223.                 carry = sival.sihigh;
  224.                 sival.ivalue = ((FULL) *h1++) * digit
  225.                     + ((FULL) *hd) + carry;
  226.                 *hd++ = sival.silow;
  227.                 carry = sival.sihigh;
  228.                 sival.ivalue = ((FULL) *h1++) * digit
  229.                     + ((FULL) *hd) + carry;
  230.                 *hd++ = sival.silow;
  231.                 carry = sival.sihigh;
  232.             }
  233.             while (len--) {
  234.                 sival.ivalue = ((FULL) *h1++) * digit
  235.                     + ((FULL) *hd) + carry;
  236.                 *hd++ = sival.silow;
  237.                 carry = sival.sihigh;
  238.             }
  239.             while (carry) {
  240.                 sival.ivalue = ((FULL) *hd) + carry;
  241.                 *hd++ = sival.silow;
  242.                 carry = sival.sihigh;
  243.             }
  244.         }
  245.  
  246.         /*
  247.          * Now return the true size of the number.
  248.          */
  249.         len = sizetotal;
  250.         hd = &ans[len - 1];
  251.         while ((*hd == 0) && (len > 1)) {
  252.             hd--;
  253.             len--;
  254.         }
  255.         return len;
  256.     }
  257.  
  258.     /*
  259.      * Need to multiply by a large number.
  260.      * Allocate temporary space for calculations, and calculate the
  261.      * value for the shift.  The shift value is 1/2 the size of the
  262.      * larger (first) number (rounded up).  The amount of temporary
  263.      * space needed is twice the size of the shift, plus one more word
  264.      * for the multiply to use.
  265.      */
  266.     shift = (size1 + 1) / 2;
  267.     temp = tempbuf;
  268.     tempbuf += (2 * shift) + 1;
  269.  
  270.     /*
  271.      * Determine the sizes and locations of all the numbers.
  272.      * The value of sizeC can be negative, and this is checked later.
  273.      * The value of sizeD is limited by the full size of the number.
  274.      */
  275.     baseA = v1 + shift;
  276.     baseB = v1;
  277.     /* 
  278.      * XXX - Saber-C Version 3.1 says:
  279.      *
  280.      *    W#26, Storing a bad pointer into auto variable dmul`baseC.
  281.      *
  282.      * This warning is issued during the regression test #026
  283.      * (read cryrand).
  284.      *
  285.      * Saver-C claims that v2+shift is past the end of allocated
  286.      * memory for v2.  When this warning is first issued, shift
  287.      * has the value 51.
  288.      *
  289.      * The call stack is:
  290.      *
  291.      *   domul(0x165ca0, 101, 0x160998, 40, 0x16d0a8) at "zmul.c":315 
  292.      *   zmul(0x2ddf88, 0x2ddf8c, 0x2ddc98) at "zmul.c":73 
  293.      *   qsqrt(0x38defc, 0x38ea94) at "qfunc.c":248 
  294.      *   qln(0x38defc, 0x38de70) at "qtrans.c":589 
  295.      *   qpower(0x38eacc, 0x38a398, 0x38a018) at "qtrans.c":657 
  296.      *   powervalue(0x1740f8,0x174118,0x195234,0x19523c) at "value.c":1009 
  297.      *   f_power(2, 0x533278) at "func.c":1188 
  298.      *   builtinfunc(117, 2, 0x174118) at "func.c":354 
  299.      *   o_call(0x5328d8, 117, 2) at "opcodes.c":2094 
  300.      *   calculate(0x5328d8, 1) at "opcodes.c":288 
  301.      *   o_usercall(0x5328d8, 54, 1) at "opcodes.c":2082 
  302.      *   calculate(0x48ffa8, 4) at "opcodes.c":288 
  303.      *   o_usercall(0x48ffa8, 53, 1) at "opcodes.c":2082 
  304.      *   calculate(0x1652a0, 1) at "opcodes.c":288 
  305.      *   o_usercall(0x1652a0, 57, 1) at "opcodes.c":2082 
  306.      *   calculate(0x16cca8, 0) at "opcodes.c":288 
  307.      *   evaluate(0) at "codegen.c":167 
  308.      *   getcommands(0) at "codegen.c":106 
  309.      *   getcommands(0) at "codegen.c":76 
  310.      *   getcommands(1) at "codegen.c":76 
  311.      *   main(-1, 0x181f8c) at "calc.c":155 
  312.      *
  313.      * Purify does not complain about this code.  Who is right?
  314.      */
  315.     baseC = v2 + shift;
  316.     baseD = v2;
  317.     baseAB = ans;
  318.     baseDC = ans + shift;
  319.     baseAC = ans + shift * 2;
  320.     baseBD = ans;
  321.  
  322.     sizeA = size1 - shift;
  323.     sizeC = size2 - shift;
  324.  
  325.     sizeB = shift;
  326.     hd = &baseB[shift - 1];
  327.     while ((*hd == 0) && (sizeB > 1)) {
  328.         hd--;
  329.         sizeB--;
  330.     }
  331.  
  332.     sizeD = shift;
  333.     if (sizeD > size2)
  334.         sizeD = size2;
  335.     hd = &baseD[sizeD - 1];
  336.     while ((*hd == 0) && (sizeD > 1)) {
  337.         hd--;
  338.         sizeD--;
  339.     }
  340.  
  341.     /*
  342.      * If the smaller number has a high half of zero, then calculate
  343.      * the result by breaking up the first number into two numbers
  344.      * and combining the results using the obvious formula:
  345.      *    (A*S+B) * D = (A*D)*S + B*D
  346.      */
  347.     if (sizeC <= 0) {
  348.         len = domul(baseB, sizeB, baseD, sizeD, ans);
  349.         hd = &ans[len];
  350.         len = sizetotal - len;
  351.         while (len--)
  352.             *hd++ = 0;
  353.  
  354.         /*
  355.          * Add the second number into the first number, shifted
  356.          * over at the correct position.
  357.          */
  358.         len = domul(baseA, sizeA, baseD, sizeD, temp);
  359.         h1 = temp;
  360.         hd = ans + shift;
  361.         carry = 0;
  362.         while (len--) {
  363.             sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  364.             *hd++ = sival.silow;
  365.             carry = sival.sihigh;
  366.         }
  367.         while (carry) {
  368.             sival.ivalue = ((FULL) *hd) + carry;
  369.             *hd++ = sival.silow;
  370.             carry = sival.sihigh;
  371.         }
  372.  
  373.         /*
  374.          * Determine the final size of the number and return it.
  375.          */
  376.         len = sizetotal;
  377.         hd = &ans[len - 1];
  378.         while ((*hd == 0) && (len > 1)) {
  379.             hd--;
  380.             len--;
  381.         }
  382.         tempbuf = temp;
  383.         return len;
  384.     }
  385.  
  386.     /*
  387.      * Now we know that the high halfs of the numbers are nonzero,
  388.      * so we can use the complete formula.
  389.      *    (A*S+B)*(C*S+D) = (S^2+S)*A*C + S*(A-B)*(D-C) + (S+1)*B*D.
  390.      * The steps are done in the following order:
  391.      *    A-B
  392.      *    D-C
  393.      *    (A-B)*(D-C)
  394.      *    S^2*A*C + B*D
  395.      *    (S^2+S)*A*C + (S+1)*B*D                (*)
  396.      *    (S^2+S)*A*C + S*(A-B)*(D-C) + (S+1)*B*D
  397.      *
  398.      * Note: step (*) above can produce a result which is larger than
  399.      * the final product will be, and this is where the extra word
  400.      * needed in the product comes from.  After the final subtraction is
  401.      * done, the result fits in the expected size.  Using the extra word
  402.      * is easier than suppressing the carries and borrows everywhere.
  403.      *
  404.      * Begin by forming the product (A-B)*(D-C) into a temporary
  405.      * location that we save until the final step.  Do each subtraction
  406.      * at positions 0 and S.  Be very careful about the relative sizes
  407.      * of the numbers since this result can be negative.  For the first
  408.      * step calculate the absolute difference of A and B into a temporary
  409.      * location at position 0 of the result.  Negate the sign if A is
  410.      * smaller than B.
  411.      */
  412.     neg = FALSE;
  413.     if (sizeA == sizeB) {
  414.         len = sizeA;
  415.         h1 = &baseA[len - 1];
  416.         h2 = &baseB[len - 1];
  417.         while ((len > 1) && (*h1 == *h2)) {
  418.             len--;
  419.             h1--;
  420.             h2--;
  421.         }
  422.     }
  423.     if ((sizeA > sizeB) || ((sizeA == sizeB) && h1 && h2 && (*h1 > *h2))) {
  424.         h1 = baseA;
  425.         h2 = baseB;
  426.         sizeAB = sizeA;
  427.         subsize = sizeB;
  428.     } else {
  429.         neg = !neg;
  430.         h1 = baseB;
  431.         h2 = baseA;
  432.         sizeAB = sizeB;
  433.         subsize = sizeA;
  434.     }
  435.     copysize = sizeAB - subsize;
  436.  
  437.     hd = baseAB;
  438.     carry = 0;
  439.     while (subsize--) {
  440.         sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
  441.         *hd++ = BASE1 - sival.silow;
  442.         carry = sival.sihigh;
  443.     }
  444.     while (copysize--) {
  445.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  446.         *hd++ = BASE1 - sival.silow;
  447.         carry = sival.sihigh;
  448.     }
  449.  
  450.     hd = &baseAB[sizeAB - 1];
  451.     while ((*hd == 0) && (sizeAB > 1)) {
  452.         hd--;
  453.         sizeAB--;
  454.     }
  455.  
  456.     /*
  457.      * This completes the calculation of abs(A-B).  For the next step
  458.      * calculate the absolute difference of D and C into a temporary
  459.      * location at position S of the result.  Negate the sign if C is
  460.      * larger than D.
  461.      */
  462.     if (sizeC == sizeD) {
  463.         len = sizeC;
  464.         h1 = &baseC[len - 1];
  465.         h2 = &baseD[len - 1];
  466.         while ((len > 1) && (*h1 == *h2)) {
  467.             len--;
  468.             h1--;
  469.             h2--;
  470.         }
  471.     }
  472.     if ((sizeC > sizeD) || ((sizeC == sizeD) && (*h1 > *h2)))
  473.     {
  474.         neg = !neg;
  475.         h1 = baseC;
  476.         h2 = baseD;
  477.         sizeDC = sizeC;
  478.         subsize = sizeD;
  479.     } else {
  480.         h1 = baseD;
  481.         h2 = baseC;
  482.         sizeDC = sizeD;
  483.         subsize = sizeC;
  484.     }
  485.     copysize = sizeDC - subsize;
  486.  
  487.     hd = baseDC;
  488.     carry = 0;
  489.     while (subsize--) {
  490.         sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
  491.         *hd++ = BASE1 - sival.silow;
  492.         carry = sival.sihigh;
  493.     }
  494.     while (copysize--) {
  495.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  496.         *hd++ = BASE1 - sival.silow;
  497.         carry = sival.sihigh;
  498.     }
  499.     hd = &baseDC[sizeDC - 1];
  500.     while ((*hd == 0) && (sizeDC > 1)) {
  501.         hd--;
  502.         sizeDC--;
  503.     }
  504.  
  505.     /*
  506.      * This completes the calculation of abs(D-C).  Now multiply
  507.      * together abs(A-B) and abs(D-C) into a temporary location,
  508.      * which is preserved until the final steps.
  509.      */
  510.     baseABDC = temp;
  511.     sizeABDC = domul(baseAB, sizeAB, baseDC, sizeDC, baseABDC);
  512.  
  513.     /*
  514.      * Now calculate B*D and A*C into one of their two final locations.
  515.      * Make sure the high order digits of the products are zeroed since
  516.      * this initializes the final result.  Be careful about this zeroing
  517.      * since the size of the high order words might be smaller than
  518.      * the shift size.  Do B*D first since the multiplies use one more
  519.      * word than the size of the product.  Also zero the final extra
  520.      * word in the result for possible carries to use.
  521.      */
  522.     len = domul(baseB, sizeB, baseD, sizeD, baseBD);
  523.     hd = &baseBD[len];
  524.     len = shift * 2 - len;
  525.     while (len--)
  526.         *hd++ = 0;
  527.  
  528.     len = domul(baseA, sizeA, baseC, sizeC, baseAC);
  529.     hd = &baseAC[len];
  530.     len = sizetotal - shift * 2 - len + 1;
  531.     while (len--)
  532.         *hd++ = 0;
  533.  
  534.     /*
  535.      * Now add in A*C and B*D into themselves at the other shifted
  536.      * position that they need.  This addition is tricky in order to
  537.      * make sure that the two additions cannot interfere with each other.
  538.      * Therefore we first add in the top half of B*D and the lower half
  539.      * of A*C.  The sources and destinations of these two additions
  540.      * overlap, and so the same answer results from the two additions,
  541.      * thus only two pointers suffice for both additions.  Keep the
  542.      * final carry from these additions for later use since we cannot
  543.      * afford to change the top half of A*C yet.
  544.      */
  545.     h1 = baseBD + shift;
  546.     h2 = baseAC;
  547.     carryACBD = 0;
  548.     len = shift;
  549.     while (len--) {
  550.         sival.ivalue = ((FULL) *h1) + ((FULL) *h2) + carryACBD;
  551.         *h1++ = sival.silow;
  552.         *h2++ = sival.silow;
  553.         carryACBD = sival.sihigh;
  554.     }
  555.  
  556.     /*
  557.      * Now add in the bottom half of B*D and the top half of A*C.
  558.      * These additions are straightforward, except that A*C should
  559.      * be done first because of possible carries from B*D, and the
  560.      * top half of A*C might not exist.  Add in one of the carries
  561.      * from the previous addition while we are at it.
  562.      */
  563.     h1 = baseAC + shift;
  564.     hd = baseAC;
  565.     carry = carryACBD;
  566.     len = sizetotal - 3 * shift;
  567.     while (len--) {
  568.         sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  569.         *hd++ = sival.silow;
  570.         carry = sival.sihigh;
  571.     }
  572.     while (carry) {
  573.         sival.ivalue = ((FULL) *hd) + carry;
  574.         *hd++ = sival.silow;
  575.         carry = sival.sihigh;
  576.     }
  577.  
  578.     h1 = baseBD;
  579.     hd = baseBD + shift;
  580.     carry = 0;
  581.     len = shift;
  582.     while (len--) {
  583.         sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  584.         *hd++ = sival.silow;
  585.         carry = sival.sihigh;
  586.     }
  587.     while (carry) {
  588.         sival.ivalue = ((FULL) *hd) + carry;
  589.         *hd++ = sival.silow;
  590.         carry = sival.sihigh;
  591.     }
  592.  
  593.     /*
  594.      * Now finally add in the other delayed carry from the
  595.      * above addition.
  596.      */
  597.     hd = baseAC + shift;
  598.     while (carryACBD) {
  599.         sival.ivalue = ((FULL) *hd) + carryACBD;
  600.         *hd++ = sival.silow;
  601.         carryACBD = sival.sihigh;
  602.     }
  603.  
  604.     /*
  605.      * Now finally add or subtract (A-B)*(D-C) into the final result at
  606.      * the correct position (S), according to whether it is positive or
  607.      * negative.  When subtracting, the answer cannot go negative.
  608.      */
  609.     h1 = baseABDC;
  610.     hd = ans + shift;
  611.     carry = 0;
  612.     len = sizeABDC;
  613.     if (neg) {
  614.         while (len--) {
  615.             sival.ivalue = BASE1 - ((FULL) *hd) +
  616.                 ((FULL) *h1++) + carry;
  617.             *hd++ = BASE1 - sival.silow;
  618.             carry = sival.sihigh;
  619.         }
  620.         while (carry) {
  621.             sival.ivalue = BASE1 - ((FULL) *hd) + carry;
  622.             *hd++ = BASE1 - sival.silow;
  623.             carry = sival.sihigh;
  624.         }
  625.     } else {
  626.         while (len--) {
  627.             sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  628.             *hd++ = sival.silow;
  629.             carry = sival.sihigh;
  630.         }
  631.         while (carry) {
  632.             sival.ivalue = ((FULL) *hd) + carry;
  633.             *hd++ = sival.silow;
  634.             carry = sival.sihigh;
  635.         }
  636.     }
  637.  
  638.     /*
  639.      * Finally determine the size of the final result and return that.
  640.      */
  641.     len = sizetotal;
  642.     hd = &ans[len - 1];
  643.     while ((*hd == 0) && (len > 1)) {
  644.         hd--;
  645.         len--;
  646.     }
  647.     tempbuf = temp;
  648.     return len;
  649. }
  650.  
  651.  
  652. /*
  653.  * Square a number by using the following formula recursively:
  654.  *    (A*S+B)^2 = (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2
  655.  * where S is a power of 2^16, and so multiplies by it are shifts,
  656.  * and A and B are the left and right halfs of the number to square.
  657.  */
  658. void
  659. zsquare(z, res)
  660.     ZVALUE z, *res;
  661. {
  662.     LEN len;
  663.  
  664.     if (ziszero(z)) {
  665.         *res = _zero_;
  666.         return;
  667.     }
  668.     if (zisunit(z)) {
  669.         *res = _one_;
  670.         return;
  671.     }
  672.  
  673.     /*
  674.      * Allocate a temporary array if necessary for the recursion to use.
  675.      * The array needs to be allocated large enough for all of the
  676.      * temporary results to fit in.  This size is about 3 times the
  677.      * size of the original number, since each recursion level uses 3/2
  678.      * of the size of its given number, and whose size is 1/2 the size
  679.      * of the previous level.  The sum of the infinite series is 3.
  680.      * Allocate some extra words for rounding up the sizes.
  681.      */
  682.     len = 3 * z.len + 32;
  683.     tempbuf = zalloctemp(len);
  684.  
  685.     res->sign = 0;
  686.     res->v = alloc((z.len+1) * 2);
  687.     res->len = dosquare(z.v, z.len, res->v);
  688. }
  689.  
  690.  
  691. /*
  692.  * Recursive routine to square a number by splitting it up into two numbers
  693.  * of half the size, and using the results of squaring the subpieces.
  694.  * The result is placed in the indicated array, which must be large
  695.  * enough for the result (size * 2).  Returns the size of the result.
  696.  * This uses a temporary array which must be 3 times as large as the
  697.  * size of the number at the top level recursive call.
  698.  */
  699. static LEN
  700. dosquare(vp, size, ans)
  701.     HALF *vp;        /* value to be squared */
  702.     LEN size;        /* length of value to square */
  703.     HALF *ans;        /* location for result */
  704. {
  705.     LEN shift;        /* amount numbers are shifted by */
  706.     LEN sizeA;        /* size of left half of number to square */
  707.     LEN sizeB;        /* size of right half of number to square */
  708.     LEN sizeAA;        /* size of square of left half */
  709.     LEN sizeBB;        /* size of square of right half */
  710.     LEN sizeAABB;        /* size of sum of squares of A and B */
  711.     LEN sizeAB;        /* size of difference of A and B */
  712.     LEN sizeABAB;        /* size of square of difference of A and B */
  713.     LEN subsize;        /* size of difference of halfs */
  714.     LEN copysize;        /* size of number left to copy */
  715.     LEN sumsize;        /* size of sum */
  716.     LEN sizetotal;        /* total size of square */
  717.     LEN len;        /* temporary length */
  718.     LEN len1;        /* another temporary length */
  719.     FULL carry;        /* carry digit for small multiplications */
  720.     FULL digit;        /* single digit multiplying by */
  721.     HALF *temp;        /* base for temporary calculations */
  722.     HALF *baseA;        /* base of left half of number */
  723.     HALF *baseB;        /* base of right half of number */
  724.     HALF *baseAA;        /* base of square of left half of number */
  725.     HALF *baseBB;        /* base of square of right half of number */
  726.     HALF *baseAABB;        /* base of sum of squares of A and B */
  727.     HALF *baseAB;        /* base of difference of A and B */
  728.     HALF *baseABAB;        /* base of square of difference of A and B */
  729.     register HALF *hd, *h1, *h2, *h3;    /* for inner loops */
  730.     SIUNION sival;        /* for addition of digits */
  731.  
  732.     /*
  733.      * First trim the number of leading zeroes.
  734.      */
  735.     hd = &vp[size - 1];
  736.     while ((*hd == 0) && (size > 1)) {
  737.         size--;
  738.         hd--;
  739.     }
  740.     sizetotal = size + size;
  741.  
  742.     /*
  743.      * If the number has only a small number of digits, then use the
  744.      * (almost) normal multiplication method.  Multiply each halfword
  745.      * only by those halfwards further on in the number.  Missed terms
  746.      * will then be the same pairs of products repeated, and the squares
  747.      * of each halfword.  The first case is handled by doubling the
  748.      * result.  The second case is handled explicitly.  The number of
  749.      * digits where the algorithm changes is settable from 2 to maxint.
  750.      */
  751.     if (size < _sq2_) {
  752.         hd = ans;
  753.         len = sizetotal;
  754.         while (len--)
  755.             *hd++ = 0;
  756.  
  757.         h2 = vp;
  758.         hd = ans + 1;
  759.         for (len = size; len--; hd += 2) {
  760.             digit = (FULL) *h2++;
  761.             if (digit == 0)
  762.                 continue;
  763.             h3 = h2;
  764.             h1 = hd;
  765.             carry = 0;
  766.             len1 = len;
  767.             while (len1 >= 4) {    /* expand the loop some */
  768.                 len1 -= 4;
  769.                 sival.ivalue = (digit * ((FULL) *h3++))
  770.                     + ((FULL) *h1) + carry;
  771.                 *h1++ = sival.silow;
  772.                 sival.ivalue = (digit * ((FULL) *h3++))
  773.                     + ((FULL) *h1) + ((FULL) sival.sihigh);
  774.                 *h1++ = sival.silow;
  775.                 sival.ivalue = (digit * ((FULL) *h3++))
  776.                     + ((FULL) *h1) + ((FULL) sival.sihigh);
  777.                 *h1++ = sival.silow;
  778.                 sival.ivalue = (digit * ((FULL) *h3++))
  779.                     + ((FULL) *h1) + ((FULL) sival.sihigh);
  780.                 *h1++ = sival.silow;
  781.                 carry = sival.sihigh;
  782.             }
  783.             while (len1--) {
  784.                 sival.ivalue = (digit * ((FULL) *h3++))
  785.                     + ((FULL) *h1) + carry;
  786.                 *h1++ = sival.silow;
  787.                 carry = sival.sihigh;
  788.             }
  789.             while (carry) {
  790.                 sival.ivalue = ((FULL) *h1) + carry;
  791.                 *h1++ = sival.silow;
  792.                 carry = sival.sihigh;
  793.             }
  794.         }
  795.  
  796.         /*
  797.          * Now double the result.
  798.          * There is no final carry to worry about because we
  799.          * handle all digits of the result which must fit.
  800.          */
  801.         carry = 0;
  802.         hd = ans;
  803.         len = sizetotal;
  804.         while (len--) {
  805.             digit = ((FULL) *hd);
  806.             sival.ivalue = digit + digit + carry;
  807.             *hd++ = sival.silow;
  808.             carry = sival.sihigh;
  809.         }
  810.  
  811.         /*
  812.          * Now add in the squares of each halfword.
  813.          */
  814.         carry = 0;
  815.         hd = ans;
  816.         h3 = vp;
  817.         len = size;
  818.         while (len--) {
  819.             digit = ((FULL) *h3++);
  820.             sival.ivalue = digit * digit + ((FULL) *hd) + carry;
  821.             *hd++ = sival.silow;
  822.             carry = sival.sihigh;
  823.             sival.ivalue = ((FULL) *hd) + carry;
  824.             *hd++ = sival.silow;
  825.             carry = sival.sihigh;
  826.         }
  827.         while (carry) {
  828.             sival.ivalue = ((FULL) *hd) + carry;
  829.             *hd++ = sival.silow;
  830.             carry = sival.sihigh;
  831.         }
  832.  
  833.         /*
  834.          * Finally return the size of the result.
  835.          */
  836.         len = sizetotal;
  837.         hd = &ans[len - 1];
  838.         while ((*hd == 0) && (len > 1)) {
  839.             len--;
  840.             hd--;
  841.         }
  842.         return len;
  843.     }
  844.  
  845.     /*
  846.      * The number to be squared is large.
  847.      * Allocate temporary space and determine the sizes and
  848.      * positions of the values to be calculated.
  849.      */
  850.     temp = tempbuf;
  851.     tempbuf += (3 * (size + 1) / 2);
  852.  
  853.     sizeA = size / 2;
  854.     sizeB = size - sizeA;
  855.     shift = sizeB;
  856.     baseA = vp + sizeB;
  857.     baseB = vp;
  858.     baseAA = &ans[shift * 2];
  859.     baseBB = ans;
  860.     baseAABB = temp;
  861.     baseAB = temp;
  862.     baseABAB = &temp[shift];
  863.  
  864.     /*
  865.      * Trim the second number of leading zeroes.
  866.      */
  867.     hd = &baseB[sizeB - 1];
  868.     while ((*hd == 0) && (sizeB > 1)) {
  869.         sizeB--;
  870.         hd--;
  871.     }
  872.  
  873.     /*
  874.      * Now to proceed to calculate the result using the formula.
  875.      *    (A*S+B)^2 = (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2.
  876.      * The steps are done in the following order:
  877.      *    S^2*A^2 + B^2
  878.      *    A^2 + B^2
  879.      *    (S^2+S)*A^2 + (S+1)*B^2
  880.      *    (A-B)^2
  881.      *    (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2.
  882.      *
  883.      * Begin by forming the squares of two the halfs concatenated
  884.      * together in the final result location.  Make sure that the
  885.      * highest words of the results are zero.
  886.      */
  887.     sizeBB = dosquare(baseB, sizeB, baseBB);
  888.     hd = &baseBB[sizeBB];
  889.     len = shift * 2 - sizeBB;
  890.     while (len--)
  891.         *hd++ = 0;
  892.  
  893.     sizeAA = dosquare(baseA, sizeA, baseAA);
  894.     hd = &baseAA[sizeAA];
  895.     len = sizetotal - shift * 2 - sizeAA;
  896.     while (len--)
  897.         *hd++ = 0;
  898.  
  899.     /*
  900.      * Sum the two squares into a temporary location.
  901.      */
  902.     if (sizeAA >= sizeBB) {
  903.         h1 = baseAA;
  904.         h2 = baseBB;
  905.         sizeAABB = sizeAA;
  906.         sumsize = sizeBB;
  907.     } else {
  908.         h1 = baseBB;
  909.         h2 = baseAA;
  910.         sizeAABB = sizeBB;
  911.         sumsize = sizeAA;
  912.     }
  913.     copysize = sizeAABB - sumsize;
  914.  
  915.     hd = baseAABB;
  916.     carry = 0;
  917.     while (sumsize--) {
  918.         sival.ivalue = ((FULL) *h1++) + ((FULL) *h2++) + carry;
  919.         *hd++ = sival.silow;
  920.         carry = sival.sihigh;
  921.     }
  922.     while (copysize--) {
  923.         sival.ivalue = ((FULL) *h1++) + carry;
  924.         *hd++ = sival.silow;
  925.         carry = sival.sihigh;
  926.     }
  927.     if (carry) {
  928.         *hd = (HALF)carry;
  929.         sizeAABB++;
  930.     }
  931.  
  932.     /*
  933.      * Add the sum back into the previously calculated squares
  934.      * shifted over to the proper location.
  935.      */
  936.     h1 = baseAABB;
  937.     hd = ans + shift;
  938.     carry = 0;
  939.     len = sizeAABB;
  940.     while (len--) {
  941.         sival.ivalue = ((FULL) *hd) + ((FULL) *h1++) + carry;
  942.         *hd++ = sival.silow;
  943.         carry = sival.sihigh;
  944.     }
  945.     while (carry) {
  946.         /*
  947.          * XXX - Saber-C Version 3.1 says:
  948.          *
  949.          *    W#113, Using allocated data what has not been set.
  950.          *
  951.          * This warning is issued during the regression test #622
  952.          * (algcheck 1).
  953.          *
  954.          * Saber-C claims that *hd has been allocated but not set.
  955.          * When this warning was first issued, hd == 0x19547c
  956.          * which is 12 bytes off from ans.  The value of carry
  957.          * was 1, the value of shift was 2, len was -1 and sizeAABB
  958.          * was 4, sizeBB as 4, sizeAA was 2.  The value stored
  959.          * in *hd was 0xbfbf.
  960.          *
  961.          * The call stack is:
  962.          *
  963.          *   dosquare(0x183e88, 3, 0x195470) at "zmul.c":981 
  964.          *   zsquare(0x1952b0, 0x4c1348) at "zmul.c":687 
  965.          *   zpowermod(0x5d1d58,0x5d1d5c,0x5d1d60,0x390218) at "zmod.c":679 
  966.          *   qpowermod(0x38a088, 0x3901a8, 0x38a964) at "qfunc.c":76 
  967.          *   builtinfunc(113, 3, 0x1740f0) at "func.c":386 
  968.          *   o_call(0x3b5de0, 113, 3) at "opcodes.c":2094 
  969.          *   calculate(0x3b5de0, 3) at "opcodes.c":288 
  970.          *   o_usercall(0x3b5de0, 13, 3) at "opcodes.c":2082 
  971.          *   calculate(0x47dfa8, 0) at "opcodes.c":288 
  972.          *   o_usercall(0x47dfa8, 14, 0) at "opcodes.c":2082 
  973.          *   calculate(0x16cca8, 0) at "opcodes.c":288 
  974.          *   evaluate(0) at "codegen.c":167 
  975.          *   getcommands(0) at "codegen.c":106 
  976.          *   getcommands(1) at "codegen.c":76 
  977.          *   main(-1, 0x181f8c) at "calc.c":155 
  978.          *
  979.          * Purify does not report this.  Who is right?
  980.          */
  981.         sival.ivalue = ((FULL) *hd) + carry;
  982.         *hd++ = sival.silow;
  983.         carry = sival.sihigh;
  984.     }
  985.  
  986.     /*
  987.      * Calculate the absolute value of the difference of the two halfs
  988.      * into a temporary location.
  989.      */
  990.     if (sizeA == sizeB) {
  991.         len = sizeA;
  992.         h1 = &baseA[len - 1];
  993.         h2 = &baseB[len - 1];
  994.         while ((len > 1) && (*h1 == *h2)) {
  995.             len--;
  996.             h1--;
  997.             h2--;
  998.         }
  999.     }
  1000.     if ((sizeA > sizeB) || ((sizeA == sizeB) && (*h1 > *h2)))
  1001.     {
  1002.         h1 = baseA;
  1003.         h2 = baseB;
  1004.         sizeAB = sizeA;
  1005.         subsize = sizeB;
  1006.     } else {
  1007.         h1 = baseB;
  1008.         h2 = baseA;
  1009.         sizeAB = sizeB;
  1010.         subsize = sizeA;
  1011.     }
  1012.     copysize = sizeAB - subsize;
  1013.  
  1014.     hd = baseAB;
  1015.     carry = 0;
  1016.     while (subsize--) {
  1017.         sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
  1018.         *hd++ = BASE1 - sival.silow;
  1019.         carry = sival.sihigh;
  1020.     }
  1021.     while (copysize--) {
  1022.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  1023.         *hd++ = BASE1 - sival.silow;
  1024.         carry = sival.sihigh;
  1025.     }
  1026.  
  1027.     hd = &baseAB[sizeAB - 1];
  1028.     while ((*hd == 0) && (sizeAB > 1)) {
  1029.         sizeAB--;
  1030.         hd--;
  1031.     }
  1032.  
  1033.     /*
  1034.      * Now square the number into another temporary location,
  1035.      * and subtract that from the final result.
  1036.      */
  1037.     sizeABAB = dosquare(baseAB, sizeAB, baseABAB);
  1038.  
  1039.     h1 = baseABAB;
  1040.     hd = ans + shift;
  1041.     carry = 0;
  1042.     while (sizeABAB--) {
  1043.         sival.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++) + carry;
  1044.         *hd++ = BASE1 - sival.silow;
  1045.         carry = sival.sihigh;
  1046.     }
  1047.     while (carry) {
  1048.         sival.ivalue = BASE1 - ((FULL) *hd) + carry;
  1049.         *hd++ = BASE1 - sival.silow;
  1050.         carry = sival.sihigh;
  1051.     }
  1052.  
  1053.     /*
  1054.      * Return the size of the result.
  1055.      */
  1056.     len = sizetotal;
  1057.     hd = &ans[len - 1];
  1058.     while ((*hd == 0) && (len > 1)) {
  1059.         len--;
  1060.         hd--;
  1061.     }
  1062.     tempbuf = temp;
  1063.     return len;
  1064. }
  1065.  
  1066.  
  1067. /*
  1068.  * Return a pointer to a buffer to be used for holding a temporary number.
  1069.  * The buffer will be at least as large as the specified number of HALFs,
  1070.  * and remains valid until the next call to this routine.  The buffer cannot
  1071.  * be freed by the caller.  There is only one temporary buffer, and so as to
  1072.  * avoid possible conflicts this is only used by the lowest level routines
  1073.  * such as divide, multiply, and square.
  1074.  */
  1075. HALF *
  1076. zalloctemp(len)
  1077.     LEN len;        /* required number of HALFs in buffer */
  1078. {
  1079.     HALF *hp;
  1080.     static LEN buflen;    /* current length of temp buffer */
  1081.     static HALF *bufptr;    /* pointer to current temp buffer */
  1082.  
  1083.     if (len <= buflen)
  1084.         return bufptr;
  1085.  
  1086.     /*
  1087.      * We need to grow the temporary buffer.
  1088.      * First free any existing buffer, and then allocate the new one.
  1089.      * While we are at it, make the new buffer bigger than necessary
  1090.      * in order to reduce the number of reallocations.
  1091.      */
  1092.     len += 100;
  1093.     if (buflen) {
  1094.         buflen = 0;
  1095.         free(bufptr);
  1096.     }
  1097.     hp = (HALF *) malloc(len * sizeof(HALF));
  1098.     if (hp == NULL)
  1099.         math_error("No memory for temp buffer");
  1100.     bufptr = hp;
  1101.     buflen = len;
  1102.     return hp;
  1103. }
  1104.  
  1105. /* END CODE */
  1106.